def codon_probs_of_parent_seq_new(
    selection_crepe, nt_sequence, branch_length, neutral_crepe=None, multihit_model=None
):
    """Calculate the predicted model probabilities of each codon at each site.

    Args:
        nt_sequence: A tuple of two strings, the heavy and light chain nucleotide
            sequences.
        branch_length: The branch length of the tree.
    Returns:
        a tuple of tensors of shape (L, 64) representing the predicted probabilities of each
        codon at each site.
    """
    if neutral_crepe is None:
        raise NotImplementedError("neutral_crepe is required.")

    if isinstance(nt_sequence, str) or len(nt_sequence) != 2:
        raise ValueError(
            "nt_sequence must be a pair of strings, with the first element being the heavy chain sequence and the second element being the light chain sequence."
        )

    aa_seqs = tuple(translate_sequences_mask_codons(nt_sequence))
    # We must mask any codons containing N's because we need neutral probs to
    # do simulation:
    mask = tuple(codon_mask_tensor_of(chain_nt_seq) for chain_nt_seq in nt_sequence)
    rates, csps = trimmed_shm_outputs_of_parent_pair(neutral_crepe, nt_sequence)

    selection_factors = selection_crepe([aa_seqs])[0]

    if selection_crepe.model.hyperparameters["output_dim"] == 1:
        # Need to upgrade single selection factor to 20 selection factors, all
        # equal except for the one for the parent sequence, which should be
        # 1 (0 in log space).
        new_selection_factors = []
        for aa_seq, old_selection_factors in zip(aa_seqs, selection_factors):
            if len(aa_seq) == 0:
                new_selection_factors.append(torch.empty(0, 20, dtype=old_selection_factors.dtype))
            else:
                parent_indices = aa_idx_tensor_of_str_ambig(aa_seq)
                # print(old_selection_factors)
                new_selection_factors.append(
                    # Selection factors are expected to be in linear space here
                    molevol.lift_to_per_aa_selection_factors(old_selection_factors, parent_indices)
                )
        selection_factors = tuple(new_selection_factors)

    parent_nt_idxs = tuple(
        nt_idx_tensor_of_str(nt_chain_seq.replace("N", "A")) for nt_chain_seq in nt_sequence
    )
    codon_probs = []
    for parent_idxs, nt_csps, nt_rates, sel_matrix in zip(parent_nt_idxs, csps, rates, selection_factors):
        if len(parent_idxs) > 0:
            nt_mut_probs = 1.0 - torch.exp(-branch_length * nt_rates)
            codon_mutsel, _ = molevol.build_codon_mutsel(
                parent_idxs.reshape(-1, 3),
                nt_mut_probs.reshape(-1, 3),
                nt_csps.reshape(-1, 3, 4),
                sel_matrix,
                multihit_model=multihit_model,
            )
            codon_probs.append(molevol.zero_stop_codon_probs(molevol.flatten_codons(clamp_probability(codon_mutsel))))
        else:
            codon_probs.append(torch.empty(0, 64, dtype=torch.float32))

    return tuple(codon_probs)